Novel approach to description of spin liquid phases in low-dimensional quantum 

ant iferromagnet s . 
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We consider quantum spin systems with dimerization, which at strong coupling have singlet ground 
states. To account for strong correlations, the S = 1 elementary excitations are described as dilute 
Bose gas with infinite on-site repulsion. This approach is applied to the two-layer Heisenberg 
antiferromagnet at T = with general couplings. Our analytic results for the triplet gap, the 
excitation spectrum and the location of the quantum critical point are in excellent agreement with 
numerical results, obtained by dimer series expansions. 
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One of the most chaUenging problems of quantum mag- 
netism is the description of transitions between phases 
with spontaneously broken symmetry and disordered 
(spin liquid) phases. The properties of the disordered 
phases are also of great interest [Q . 

A variety of quantum spin models have been intro- 
duced in connection with the high- Tc cuprates and other 
recently discovered compounds. Examples include the 
Heisenberg ladder , the two-layer Heisenberg model ||] 
and 2D square lattice models with dimerization Q . In all 
of the above the Hamiltonian favors singlet formation of 
the spins between the chains (layers) or on neighboring 
sites. For this class of models the disordered phase is rel- 
atively well understood, since the lowest excitation above 
the singlet is a massive triplet. Another example is the 
CaV^Og lattice pi, where the spins form a singlet state 
on a plaquette. There also have been suggestions that 
dimerization of different kinds may occur in the Ji — J2 
model |,|. 

All of the models mentioned above, except for the lad- 
der, exhibit a quantum phase transition from a disor- 
dered dimer phase to a coUinear Neel phase with long 
range order in the ground state as the dimerization de- 
creases. This transition occurs due to competition be- 
tween singlet formation and antiferromagnetic order. A 
useful approach to the description of the disordered phase 
is the bond operator representation for spins, introduced 
by Chubukov g and Sachdev and Bhatt |^. This rep- 
resentation can be considered as the analog of the usual 
Holstein-Primakoff transformation for phases with un- 
broken spin rotational symmetry. Let us consider two 
S = 1/2 spins Si, S2 and introduce operators for cre- 
ation of a singlet s^O >= A={\ ti> — | iT>) and three 

triplet states tjj,a = x,y,z above a fictitious vaccum 

|o >: 4|o >= -^(1 TT> -I U>),4lo >= *7f(l TT> 

+ 1 ii>),4|0 >= ^(1 n> +1 iT>)- Then the following 
representation is exact 0: 

^"2 = ^(is^ic ± tis - iec^pjtl^) (1) 

The four operators satisfy the usual bosonic commuta- 
tion relations. In order to ensure that the physical states 



are either singlets or triplets one has to impose the con- 
dition: s^ s + t\,ta = 1. For a lattice spin system, the 
constraint is typically taken into account in a mean-field 
fashion, i.e. it is not strictly satisfied on every site, but 
only on average Q. A slightly different representation 
can be obtained by choosing the singlet as the ground 
state. Then Eq.(l) is still valid, but the operator s has 

the form: s = y 1 — tata, which formally is the resolu- 
tion of the constraint [plpl . Again, the form of s ensures 
that only physical states are present. However, it is very 
difficult to take the s term into account due to its non- 
linear nature. Expansions of the square root to infinite 
order have been proposed |g| . Unfortunately, there is no 
small parameter in this expansion and therefore the sum- 
mation is ambiguous and technically complicated. Alter- 
natively, one can use numerical techniques, based on the 
Gutzwiller projection method ||9[]. 

In this Letter we present an effective analytical method 
to deal with the hard core constraint. This approach can 
be applied to any model, for which the excitations in 
the disordered phase are triplets above a strong coupling 
singlet ground state. For definiteness we consider the 
model: 

H = J 2_^ Sii.Sij + XJ 2_^ S2i-S2j + J^y ^ Sii.S2i- 



<i,3> 



<i,3> 



(2) 



All the spins arc 1/2 and the couplings are antiferromag- 
netic (J, J_L > 0). The spins Su, 821 represent two planes 
of Heisenberg spins, coupled through the third term in 
(2). The summation in each plane is over nearest neigh- 
bors on a square lattice. In the present work we consider 
two cases: A = 1, which corresponds to the two-layer 
Heisenberg model, and A = 0, describing free spins in 
one of the planes. The latter model is interesting be- 
cause of its connection to the Kondo lattice model (at 
half filling) with an additional repulsive Hubbard inter- 
action between the conduction electrons. In the limit 
when the repulsion is strong, the charge degrees of free- 
dom are frozen, while the spin part is described by the 
Heisenberg Hamiltonian, leading to (2) (at A = 0). A 



simplified version of this model was introduced by Do- 
niach [|l^ to study the competition between local singlet 
formation (the Kondo effect) and the induced magnetic 
(RKKY) interaction between the free spins. The early 
mean- field treatment of Doniach predicted a critical J±^ , 
below which the free spins order antiferromagnetically. 
Recently this model was also studied numerically by Mat- 
sushita, Gelfand and Ishii ll^] who also found a finite 
transition point. 

For J± >> J interplane singlets are favored and the 
wave function is a product of on-site dimers. The ex- 
citations above this strong coupling ground state are 
triplets. In order to obtain the effective Hamiltonian for 
the triplets we pair the spins into inter-plane singlets by 
using (1). Alternatively, instead of applying the trans- 
formation (1), one could use perturbation theory in the 
"hopping" J, and calculate matrix elements of the type: 

1/4, < tai,ti3j\Sii.Sij\tyi,tsj >= l/4:{Sa5SyP — 6af}S^s), 

etc. The latter method is more useful when additional 
degrees of freedom are present in the problem, e.g. holes. 
For a start we neglect the constraint completely (i.e. for- 
mally set s = 1 in (1)) and obtain the effective Hamilto- 
nian: 



H — H2 + H3 + H4, 



H2 — 2_^ ^k^ka^ka 



k,a 



Sk 

2 



C^lko+h-c.) 



H.= 



(A-l)J 



<i,j>,a(3^ 



{[«ea/37iL^]3j*7i+ll-' 



(3) 

(4) 

(5) 



Ha 



(1 + A)J 



X! \}'^ait\jifii'tao " tlitljiflttpj j 



<JJ>,Q/3 



(6) 



The coefficients in (4) are: Ak = J± + (1 + A) J^k, -Sk = 
(1 -I- A) J^k, where ^k ~ {cos{kx) + cos{ky))/2. By using 
the Bogoliubov transformation tka = Wk^ka + ^ki_ka ^^ 
obtain for the excitation spectrum at the quadratic level 
{H2 only): tj^ = A^ — B^. The gap A = LUn,^^ is non- 
zero for J± > {J±)c — 2(1 + A) J and vanishes at (J_l)c, 
signaling a transition to a Neel ordered phase. The lo- 
cation of the critical point at this level of approximation 
(Jj_)^ = 2 J (A = 0), 4J (A = 1) differs significantly from 
the recent numerical results 1.39 ( |ll|, this work) and 
2.54 p^ , respectively. Let us mention that spin wave 
theory in the ordered phase works rather poorly for this 
problem, predicting for A = 1, (J^)c ~ 4.3 J [^. 

We find, in agreement with previous work MM, that 
the effect of the terms H3 and H4 on the spectrum is 



quite small and therefore can not explain the numerical 
results. We treat these terms later perturbatively. 

The dominant contribution to the renormalization of 
the spectrum comes from the constraint that only one of 
the triplet states can be excited on every site: iJ^j^L = 
0. This hard-core condition can be taken into account 
by introducing an infinite on-site repulsion between the 
bosons: 

Hu^uJ2tLthtp^t»^, U^oo (7) 

z,a/3 

Since the interaction is infinite, one has to find the ex- 
act scattering amplitude for the triplets. Our treatment 
is similar to the one used for Fermi gas with hard core, 
which appears in the theory of nuclear matter and ^He. 
The approach was initiated by Brueckner [Of . The scat- 
tering vertex Tap^.ysi'K.), K = (k, w) in the ladder ap- 
proximation satisfies the Bethe-Salpeter equation, shown 
in Fig. la. It depends on the total energy and momen- 
tum of the incoming particles K — Ki -I- K2 and has the 
structure Tap^-yS = ^{Sa-fSiss + SasSp-f). Since the inter- 
action is local and non-retarded, the equation for F can 
be readily solved with the result: 



F(K) 




g(q)g(k-q; 



^q"Lq 



-UJa 



^k-c 



(8) 



Here G(Q) is the normal triplet Green's function, i.e. 
G(k,t) ^ -i < r(tkQ(04a(0)) > and the Bogoliubov 
coefficients u^,v^ — ±1/2 -I- Ak/2tJk- The imaginary 
part of F is determined by the rule uj —f lu + iS. 

The basic approximation made in the derivation of 
r(K) is that we neglect all anomalous scattering ver- 
tices, which are present in the theory due to the exis- 
tence of anomalous Green's functions, Ga(k, f) = — i < 
'^i^-kai^)^i.ai^)) >■ We have also derived the complete 
set of equations by taking all vertices into account. How- 
ever, our key observation is that all anomalous contribu- 
tions are suppressed by an additional small parameter, 
present in the theory - the density of triplet excitations 
"^ = Ea < 4^- >= 3iV-i ^q «2 » 0.1 at J^/J » 2.5. 
We find that rii is quite small throughout the disordered 
phase, even close to the transition point. Thus the triplet 
excitations behave as a dilute, strongly interacting Bose 
gas. Consequently, since an insertion of an anomalous 
Green's function into the intermediate states of the lad- 
der in Fig. la brings powers of Wq into the equation for the 
amplitude, its contribution is small. Therefore Eq.(8) can 
be considered as the first term in an expansion in powers 
of the gas parameter n^. To be consistent, we also neglect 
the second term in (8), since it contains Uq. 

The self-energy, corresponding to the scattering am- 
plitude r is found as a sum of the diagrams shown in 
Fig.lb: 



q 



q,UJ-Uq). 



(9) 



Let us stress again that at our level of approximation 
(dilute gas), there is only a normal self-energy. Next, 
in order to find the renomalized spectrum, one has to 
solve the coupled Dyson equations for the normal and 
anomalous Green's functions. Since the procedure is well 
know from the theory of a Bose gas, we only write the 
final result for the normal Green's function mSi' 



GiKu) 



■ A]^ + T,{k, -u) 



[uj + Ai, + I](k, -uj)] [uj-Ai,- S(k, Lj)] + Bl 

(10) 



After separating this equation into a quasiparticle con- 
tribution and incoherent background, we find: 



G(k,c.)= ^•'^^ 



^'^'^% + G,„,. (11) 



u) — ilk + J(5 Lo -\- flk — iS 

The renormalized triplet spectrum and the renomaliza- 
tion constant are: 



fik = ^kx/(Ak + S(k,0))2 



9S\ 



Bl 



(12) 



(13) 



The renomalized Bogoliubov coefficients in (11) are: 



Ui,V^ = ±^ 



1 , Zk(Ak + S(k,0)) 



2f7i< 



(14) 



Equations (8,9,12-14) have to be solved self-consistently 
for 5](k, 0)) and Zk. From Eq.(ll) it is also clear that 
one has to replace Uk — > ^/Z^Uu, Wk -^ V^kVic in (8) and 
(9) (and also in (15,16), see below). 

We have found that effect of H^ and H4 on the quasi- 
particle spectrum is small, compared to the renormal- 
ization due to Hij. However, these two terms have to 
be included for the precise determination of the critical 
point. We treat -ff4 in mean field theory, by splitting the 
quartic operator products into all possible pairs. This is 
equivalent to taking only one-loop diagrams (first order 
in J) into account. These diagrams renormalize the two 
coefficients: 



^k ^ ^k + (1 + A) JCk 



nT.^^<^ 



i?k^i?k-(i + A)./a^^eq«q 



q 



(15) 



(16) 



This concludes the solution of the two-layer problem (A 



To solve the case A = we also have to take into 
account H3. It is convenient to rewrite H3 in terms of 
the Bogoliubov transformed operators fka , ^^ ' since in 
this way only the normal Green's functions remain. To 
one loop order ( J^) the renormalization of the spectrum 
is determined by the sum of the two diagrams in Fig.lc. 
The formula for the interaction vertex in Fig.lc. is quite 
lengthy and we do not present it here. Once the vertex is 
known, the self-energy of Fig.lc. can be easily computed, 
leading to renormalization of A^ and B^. 

The results of the self-consistent numerical solution 
are summarized in Figures 2 and 3. Figure 2 shows the 
triplet gap A = rijr.-rr a-s a function of the interlayer cou- 
pling. The transition into the Neel ordered phase occurs 
at {J±/J)c = 2.57(A = 1),1.37(A = 0). We have also 
calculated the gap by using dimer series expansions |lj] 
up to order 11(10) for A — 1(A = 0). The critical points 
are found at: {J±/J)c = 2.52(2)(A = 1), 1.39(4)(A = 0), 
or 2.537(5)(A = 1), 1.393(8)(A = 0) by fixing the critical 
exponent v = 0.71 |1^. The agreement between the an- 
alytic method and the dimer series results is excellent. 
Such a good agreement is better than might have been 
expected. Our analytic method involves approximations 
and an error of a few percent is always expected. The gap 
critical exponent v, defined as: A ^^ {J — JcY ^s v = 0.5 
in our analytical calculation, while the dimer series gives 
V ~ 0.7, in agreement with the 0(3) non-linear sigma 
model prediction. Recall that the mean field approxi- 
mation gives v — 1. Our diagrammatic approach is not 
valid very close to the critical point since the neglected 
terms in r(K) are of the form ^ ''^q/wq and thus loga- 
rithmically diverge at criticality. However this happens 
only very near to the critical point. 

The comparison of the excitation spectra, presented 
in Figure 3, shows that the agreement is very good over 
almost the whole Brillouin zone. The disagreement be- 
tween the two curves is largest at k = 0, where it is about 
5%. 

In conclusion, we have presented an effective analytical 
approach to take into account the hard core constraint 
which appears in the bond operator description of the 
dimer phase. The triplet excitations are described as a di- 
lute Bose gas with infinite on-site repulsion. We find that 
the spectrum is renormalized mostly due to the hard core, 
while the additional three and four point interactions are 
comparatively weak and can be treated perturbatively. 
The advantages of our formulation are that it is simple 
and captures the essential physics, being in agreement 
within a few percent with results obtained by dimer series 
expansions. Obvious other applications of the method 
include the 2D Heisenberg model with dimcrization, the 
Heisenberg ladder and the Kondo lattice model Il3] . The 
method can also be easily generalized to describe phases 
with spontaneously broken symmetries and nonzero tem- 
perature. 

We wish to thank M. Kuchiev and R. Eder for stimulat- 



ing discussions. The financial support of the Australian 
Research Council is gratefully acknowledged. 



FIG. 1. (a) Equation for the scattering amplitude F. (b) 
Diagrams for the self-energy, corresponding to F. (c) One 
loop diagrams, arising from the three-point interaction. 
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FIG. 2. Triplet gap as a function of interlayer coupling for 
A = (left curves) and A = 1 (right curves). The dashed 
fines with the solid circles are the results of the self-consistent 
solution. Open squares (with error bars) are from direct Pade 
approximants to the dimer series while solid lines are from 
approximants which assume v — 0.71. 



FIG. 3. Triplet excitation spectrum for A = 1 along high 
symmetry directions in the Brillouin zone. The dashed line 
with solid circles is the self-consistent solution while the solid 
lines are from direct summation of the dimer series. The 
upper (at k = 0) dimer series curve corresponds to the critical 
(within the error bar) spectrum {J±/J = 2.54), while the 
upper analytical curve is computed at J±/J = 2.6, in order 
to have the same gap A/J± = 0.05. The lower (at k = 0) 
curves correspond to J±/J = 3.33. 
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